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Abstract 

The state of a 2-D random resistor network, resulting from the simultaneous 
evolutions of two competing biased percolations, is studied in a wide range of 
bias values. Monte Carlo simulations show that when the external current I 
is below the threshold value for electrical breakdown, the network reaches a 
steady state with a nonlinear current-voltage characteristic. The properties of 
this nonlinear regime are investigated as a function of different model param- 
eters. A scaling relation is found between < R > / < R >o and I/Io, where 
< R > is the average resistance, < R >o the linear regime resistance and 
Iq the threshold value for the onset of nonlinearity. The scaling exponent is 
found to be independent of the model parameters. A similar scaling behavior 
is also found for the relative variance of resistance fluctuations. These results 
compare well with resistance measurements in composite materials performed 
in the Joule regime up to breakdown. 
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I. INTRODUCTION AND MODEL 



Electrical breakdown of disordered media has been widely studied in the last twenty years 



[M- [10]. This is due to its relevant implications on material technology and on fundamental 
aspects related to the understanding of the response properties of disordered systems to high 
external stresses PI0]. It is well known that the application of a finite stress (electrical or 
mechanical) to a disordered material generally implies a nonlinear response, which leads to 
a catastrophic behavior in the high stress limit [|I|,^] . Percolation theory provides a powerful 
approach for studying breakdown phenomena of disordered media [TT,|T^[. In this framework, 
several models have been proposed to describe the electrical breakdown of granular metals 
and of conductor-insulator composites in terms of critical phenomena near the percolation 
threshold @]- 0. Furthermore, the critical exponents characterizing these behaviors have 
been both theoretically calculated and measured in severals materials [[§]- @, ||T0 



Nevertheless, few attempts have been made so far to describe the behavior of disordered 
media over the full range of applied stress The understanding of breakdown phenomena 
in the full dynamical regime is thus unsatisfactory to the present time. On the other hand, 
important information is expected from such a study, like precursor effects, role of disorder, 
predictability of breakdown, etc. |7||J . 

The aim of this paper is to present a percolative model of sufficient generality to address 
the above issues. To this purpose, we consider a random resistor network (RRN) [|ll[ in 
which two competing percolations are present, defect generation and defect recovery, which 
determine the values of the elementary network resistances. Both processes are driven by 
an external current and by the heat exchange between the network and the thermal bath. 
Monte Carlo simulations are performed to explore the network evolution in a wide range 
of bias values. A stationary state or an irreversible breakdown of the RRN can be reached 
depending on the value of the applied current. By focusing on the steady state, we study the 
resistance and the resistance noise properties. We found that the average network resistance 
and the relative resistance noise scale with the ratio of the applied current to the current 
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value corresponding to the onset of nonlinearity. These results are discussed in connection 
with measurements in composite materials and in conducting polymers P,|T3||. 

We consider a two-dimensional, square-lattice RRN of total resistance R, made of N to t 
resistors with resistance r n . We take a square geometry N x N, where N determines the 
linear size of the network. A constant current / is applied through electrical contacts realized 
by perfectly conducting bars at the left and right sides of the network. As a consequence, 
a current i n is flowing in the nth resistor. The RRN interacts with a thermal bath at 
temperature To and the resistances r n are taken to depend linearly on the local temperature 
T n , according to: 

r n (T n ) = r [l + a(T n - T )} (1) 

where r is the resistance value of the elementary resistor at the temperature T and a is 
the temperature coefficient of the resistance. The local temperatures are calculated as in 



Ref. [10 



T n = T + A[r n tl + — — M - »•»£)]■ ( 2 ) 

In this expression, N ne i g is the number of first neighbours of the nth resistor, the parameter 
A, measured in (K/W), describes the heat coupling of each resistor with the thermal bath 
and it determines the importance of Joule heating effects. The parameter B is taken to be 
equal to 3/4 to provide a uniform heating in the perfect network configuration. We notice 
that Eq. (Q) implies an instantaneous thermalization of each resistor at the value T n , and 
then, by adopting Eq. @, we are neglecting for simplicity time dependent effects which are 
discussed in Ref. |5[]. 

In the initial state of the network, I = 0, T n = T and all the resistors are identical 
r n = r . The evolution of the RRN arises from the presence of two competing percolations. 
The first consists of generating fully insulating defects (broken resistors). This process occurs 
with probability: 

W D = exp[-E D /K B T n ] (3) 



where Ed is a characteristic activation energy and Kb the Boltzmann constant fllQ| . The 
second percolation consists of recovering the insulating defects. This process occurs with 
a probability Wr expressed as in Eq. (Q) but with a different activation energy Er. As 
a result of the competition between these two percolations, depending on the parameters 
which specify the physical properties of the system and depending on the external conditions 
(bias current and bath temperature), the RRN can reach a steady state or the percolation 
threshold. In the first case, R fluctuates around its average value < R >, while in the 
second case, an irreversible breakdown occurs, i.e. R diverges due to the existence of at 
least one continuous path of defects between the upper and lower sides of the network 
||1 1|| . By focusing on the effect of the external current, we define as the greatest current 
value for which the RNN is stationary. We notice that for biased percolation the following 
expression, E R < E D + A; B Tln[l + exp(—E D /K B T )}, represents a necessary condition for 
the existence of a steady state |L4],[l5|]. Monte Carlo simulations are performed according to 
the following procedure: (i) starting from the perfect lattice with given local currents and 
temperatures, i n and T n , (ii) resistances r n are changed according to Eq. (|TJ) and defects are 
generated with probability Wd] (in) the currents i n are calculated by solving Kirchhoff's 
loop equations; the local temperatures are updated; (iv) the temperature dependence of 
the resistances r n is again accounted for and defects are recovered with probability Wr; (v) 
R, i n and T n are finally calculated and this procedure is iterated from (ii) until electrical 
breakdown or steady state is achieved. In the last case the iteration runs long enough to 
allow a fluctuation analysis to be performed. Each iteration step can be associated with an 
elementary time step on an appropriate time scale (to be calibrated with experiments). As 
reasonable values of the parameters, simulations have been performed by taking: N = 75 
(except when differently specified), T = 300 (K), a = 10~ 3 (K" 1 ), A = 5 x 10 5 (K/W), 
E D = 0.17 (eV). Several values of E R and r have been considered: 0.026 < E R < 0.16 (eV) 
and 1 < tq < 10 (Q). The values of the external current range between 0.001 and 3 (A). 
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II. RESULTS 



We show in Fig. 1 a picture of the RRN near the electrical breakdown. In this case we 
have taken: N = 45, r = 1 (ft), E R = 0.10 (eV) and / = 0.5 (A), i.e. I > I b = 0.45 (A). 
The different levels of gray correspond to different values of r n . We can clearly see that, 
with respect to the initial state (perfect network), the network has evolved to a disordered 
state associated with the formation and growth of a channel of broken resistors elongated 
in the direction perpendicular to the applied current. This kind of damage pattern well 
reproduces the experimental pattern observed in metallic lines failed as a consequence of 
electromigration ||16|| . Typical evolutions of R are shown in Fig. 2 at increasing bias values. 
In this case N = 75 while all the other parameters are the same of Fig. 1. Thinner curves 
refer to steady state regime while the thicker curve refers to a RRN undergoing electrical 
breakdown (I > If, = 0.75). In the steady state two regimes can be identified: an Ohmic 
regime (lower two curves) and a nonlinear regime characterized by a significant increase 
of resistance (remaining curves). By focusing on the steady state regime, we report in 
Fig. 3 the average resistance < R > as a function of the applied current. The different 
curves correspond to different values of ro, i.e. to RRNs of different initial resistance, while 
the recovery activation energy is Er = 0.026 (eV). Each value have been calculated by 
considering the time average on a single realization and then making the ensemble average 
over 20 realizations. At low biases the average resistance takes a constant value < R >o 
which represents the intrinsic linear response property of the network (Ohmic regime). When 
/ is above a threshold value Jo (threshold for the nonlinearity onset), < R > increases with 
bias until the applied current reaches the I b value, above which the RRN undergoes electrical 
breakdown. Thus, in the following, we indicate with < R >f, the average value of R at / = It,, 
i.e. the last stable value of the resistance. Figure 3 also shows that by increasing ro and thus 
the initial network resistance, both J and h decrease. Precisely, we have found: I b ~ R~ s 
and Iq ~ Rq S with 5 = 0.51 ± 0.01. Therefore the ratio h/Ia = 28 ± 1 is independent of the 
initial network resistance. Moreover, we have also found < R >& / < R > = 1.85 ± 0.08. 
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The effect of the recovery activation energy on the steady state is shown in Fig. 4, which 
reports the ratio <i?>/<i?>oasa function of the applied current for different values 
of Er (in this case all the curves are obtained for ro — 1 (ft)). The overall behavior is 
similar to that shown in Fig. 3: an Ohmic regime at low bias is followed by a nonlinear 
regime for I > I . Moreover, by increasing Er both I and If, decrease and the stability 
region is thus strongly reduced. Nevertheless, an important difference between the effect of 
varying the initial network resistance and that of varying Er is that, in the last case, the 
ratio h/Io, exhibits a significant dependence on Er, as shown in Fig. 5. To investigate the 
existence of scaling relations and their universality, Fig. 6 reports the log-log plot of the 
relative variation of the average resistance, (< R > — < R >o)/ < R >o> as a function of 
the ratio I/Iq for different values of r and Er. The plot shows that all these curves collapse 
onto a single one and that the relative variation of < R > as a function of I/Iq exhibits a 
power law behavior. We conclude that, the average resistance follows the scaling relation: 

Tj£r = 9tl/h)> g(I/I )^l + (I/Io) 9 (4) 

with the scaling exponent 9 = 2.1 ± 0.1 being independent of both the initial resistance of 
the RRN and the recovery activation energy. Other simulations, performed on rectangular 
networks, confirm the same value for 9. The quadratic dependence of < R > on J, found here, 
can be understood in the spirit of mean-field theory when we consider that AR fa oRqAT 
and AT oc AR I 2 . Moreover, recalling previous results concerning the ratio I/Iq, Eq.([|) 
explains the independence of the ratio < R >& / < R >q on the initial RRN resistance (Fig. 
3) and, by contrast, its significant dependence on Er, as shown in Fig. 4. All these results 
well agree with recent measurements in the Joule regime of carbon high-density polyethylene 
composites reported in Ref. ||. 

The resistance fluctuations are analyzed for different values of Er and r . Figure 7 
reports the relative variance of resistance fluctuations, S =< AR 2 > / < R ^> , as a 
function of the external current. Curves 1, 2 and 3 (with full circles) show S for r = 1 (fl) 
and Er = 0.060, 0.043, 0.026 (eV) respectively, while the curves belonging to the set 3 are 
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obtained for Er = 0.026 (eV) and different values of r . Figure 7 points out the existence 
of two different noise regimes. The first regime occurs for / < J , i.e. when Joule heating 
effects are negligible. This noise arises from two random percolations and represents an 
intrinsic noise of the RRN, depending only on the values of Ed and Er |TJj]. The second 
regime occurs when / > I and the value of £ is found to become strongly dependent on the 
external current. By plotting £/£o as a function of I/Iq we have found that all the data of 
Fig. 7 collapse onto a single curve, as shown in Fig. 8. Moreover, a power law behavior is 
observed in the pre-breakdown region. Therefore, we can conclude that the relative variance 
of resistance fluctuations follows the scaling relation: 

£ = /(///„), /(///„)* 1 + (W (5) 

where the scaling exponent is r\ = 4.1 ± 0.1. This value of r] agrees with the values obtained 
from electrical noise measurements in conducting polymers |13 |. 

In conclusion, we have studied by Monte Carlo simulations the stationary regime of 
RRNs resulting from the simultaneous evolutions of two competing percolations. The two 
percolations consist of generating (recovering) fully insulating defects which are driven by 
an external current and by the heat exchange with a thermal bath. We have analyzed the 
behavior of the average resistance and of the relative variance of resistance fluctuations over 
a wide range of the applied current and as a function of different model parameters. We have 
found that both these quantities follow a scaling relation in terms of the ratio between the 
applied current and the current value corresponding to the nonlinearity onset. Both scaling 
exponents are found to be independent of the model parameters. These results compare 
well with resistance measurements in composite materials performed in the Joule regime up 
to breakdown II and with noise measurements in conducting polymers [ 13 . 
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FIGURE CAPTIONS 



Fig. 1 - Network evolution near the electrical breakdown. Different gray levels from black 
to white correspond to increasing values of r n from 1 to 3 (Q). 

Fig. 2 - Resistance as a function of time at increasing bias values. Thinner curves correspond 
to steady state regime and they are obtained, going from bottom to top, for I = 0.01, 0.05, 

0. 10, 0.35, 0.70, 0.75 (A). The thicker curve is obtained for I = 0.78 (A) and it corresponds 
to a RRN undergoing electrical breakdown. 

Fig. 3 - Normalized average resistance versus external bias. We take E R = 0.026 (eV) and 
E D = 0.167 (eV), while the value of r ranges between 1 Q and 10 Q. 

Fig. 4 - Normalized average resistance versus external bias. We take r = 1 (fi), E D = 0.167 
(eV) while the values of E R range between 0.026 (eV) and 0.155 (eV). 

Fig. 5 - Plot of the ratio h/h as function of the recovery activation energy. The curve is a 
fit with a power law: h/Io ~ Er°' 36 - 

Fig. 6 - Log- log plot of the relative variation of resistance versus I/Iq. Data shown in this 
figure are the same of those reported in Figs. 3 and 4. 

Fig. 7 - Relative variance of resistance fluctuations as a function of the external bias. Curves 

1, 2, 3 refer to r = 1.0 (D) and E R = 0.060, 0.043, 0.026 eV, respectively. The five curves 
belonging to set 3 are obtained with r = 1.0 (Q) (full circles), 2.5 (Q) (open squares), 5.0 
(Q) (full triangles), 7.5 (Q) (open triangles), 10.0 (Q) (full diamonds). 

Fig. 8 - Log-log plot of the relative variance of resistance fluctuations normalized to the 
same quantity calculated in the linear regime versus I / Iq. A power law fit is shown in the 
pre-breakdown regime. 
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